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Summary. — Although recent theoretical and experimental progress have consid- 
erably clarified pairing mechanisms in spin 1/2 fermionic superfluid with equally 
populated internal states, many open questions remain when the two spin popula- 
tions are mismatched. We show here that, taking advantage of the universal behavior 
characterizing the regime of infinite scattering length, the macroscopic properties of 
these systems can be simply and quantitatively understood in the regime of strong 
interactions. 



1. — Introduction 

Pairing lies at the core of the standard Bardeen-Cooper-SchriefTer mechanism for 
metal superconductivity, and the very natural question to know whether it could survive 
population imbalances between the two spin states naturally arose very soon after its 
development [TJ [2] . It was pointed out that pairing was indeed robust to some amount 
of mismatch between the chemical potentials of the two species, but the fate of the 
system after the critical imbalance is reached has long been a mystery. The absence of 
clear answer to this problem was due in particular to the absence of an experimental 
system on which the various scenarios envisioned could be tested: existence of a spatially 
modulated order parameter (Fulde, Ferrel, Larkin and Ovshinikov, or FFLO, phases) 
[3j [H [5], or the extension to trapped systems [71 [8j [9], deformed Fermi surfaces [TO] , 
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interior gap superfluidity [IT] , phase separation between a normal and a superfluid state 
through a first order phase transition p~2l [HI [15], BCS quasi-particle interactions 
[16J or onset of p-wave pairing [17] • When the strength of the interactions is varied, a 
complicated phase diagram mixing several of these scenarios is expected [18j QjJJ [20] . 

This issue was revived by the possibility of obtaining fermionic superfluids with ultra 
cold atoms [2TJ [221 [231 [211 123 12S] > where spin imbalance could be controlled and main- 
tained for a long time. This led to a series of experiments performed at MIT [271 EE] 
and Rice [29 1 30J which clearly demonstrated a phase separation between regions char- 
acterized by different polarizations (i.e. spin population imbalances, by analogy with 
magnetism). The number of phases obtained by the two groups is however different. In 
Rice experiment, the cloud is constituted of a core where both spin populations are equal, 
surrounded by a shell of majority atoms only while at MIT a third phase mixing both 
species with unequal densities is intercalated between the previous ones, a discrepancy 
which is not yet fully explained [H1[33[33[33[M 

In what follows we wish to explore the various consequences of these experiments. By 
contrast to most recent works on the subject, we would like to avoid the use of BCS mean 
field, which is known to give good qualitative insight to the problem under study, but fails 
when precise quantitative estimates are needed. Our scheme is based on a combination 
of exact variational analysis and Monte-carlo simulations. We will demonstrate that, in 
agreement with MIT experiments, three phases are expected in homogeneous systems. 
To compare with experimental results, we will make use of Local Density Approximation 
(LDA) which leads to quantitative agreement with MIT's data. Finally, following [31], we 
will show how Rice's apparently contradictory results can be interpreted as a breakdown 
of local density approximation in elongated traps. 

2. — Universal phase diagram of a homogeneous system 

Let us first consider an ensemble of spin 1/2 fermions of mass m trapped in a box of 
volume V. In the s-wave approximation, the hamiltonian H is given by 



where e& = h 2 k 2 /2m, ak,a annihilates a particle of spin a and momentum k and 
is the coupling constant characterizing s-wave interactions between atoms. This choice 
of interaction potential is singular and yields unphysical results and to get rid of the 
divergencies resulting by the zero range of the potential, we introduce an ultraviolet cut- 
off q c in momentum space (or equivalently, we work on a lattice of step l/q c )- When 
q c goes to infinity, the Lippmann-Schwinger formula obtained by the resolution of the 
two-body problem yields the following relationship between the bare coupling constant 
and the scattering length a 
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(2 ) I = ^_IyI 

V } 9b 4^ 2 a V ^e fc ' 

where the sum over k is restricted to k < q c . 

To anticipate the analysis of inhomogeneous systems, we work in the grand canonical 
ensemble, where the atom numbers fluctuate and only their expectation values are kept 
constant. Introducing the chemical potentials fi^i as Lagrange multipliers associated 
with the constraints on atom numbers, we need to find the ground state of the grand 
potential ft given by 



(3) 5 = 3-^-^. 

In what follows, we replace the minimization condition on Q = (O) by a maximization 
problem on the pressure P, using the thermodynamical relation Q = —PV. Moreover, 
we assume ^ > ^ and we restrict ourselves to the unitary limit where a = oo. This 
choice of scattering length leads to a deep simplification of the formalism, due to the 
universality characterizing this regime. Indeed, from dimensional analysis [39], we can 
show that for an arbitrary scattering length, the pressure P of a given phase is given by 
some relation 



P(ra, ft, a, /i T , fi x ) = PoOt> K ™)f{iii/ii h l//c Ft a), 

where Po is the pressure of an ideal Fermi gas with chemical potential ^ and kp] is 
the Fermi wave vector associated with /if. At unitarity, 1/k^a = and / is therefore 
function oirj = jij ji^ yielding the universal relation 



(4) TT = 9ijD, 

^0 

where #(/V/i T ) = fivi/^h )- 

Although the general minimization of the grand potential is an extremely challenging 
and still open problem, we first note that two exact eigenstates of the system can be 
found. 

1. Fully polarized ideal gas. If we consider a fully polarized system containing no 
minority atom, the interaction term in H disappears, and we are left with a pure 
ideal gas of majority atoms. The pressure of this normal phase is simply the Fermi 
pressure, and we have in particular P/Po = 1. 

2. Fully paired superfluid. Let |SF) M be the ground state of the balanced potential 
Q f = H — fi(N^ -\-N±). Since Q f commutes with the atom number operators, |SF) M 
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can be chosen as an eigenstate of both iVj^, with jV||SF) M = iVjSF)^. Going back 
to the unbalanced problem, we write ft as 



(5) Q = H + ^p^t + »l) + ^V^ T " % 

We see readily that for /i = (/xi +/i2)/2 we have ^|SF) M = Q'lSF)^, which proves 
that |SF) M is also an eigenstate of the imbalanced grand potential. The pressure 
in this superfluid phase can be calculated using known results for the unitary bal- 
anced superfluid for which the universal relationship between chemical potential 
and density reads 

(6) MT = Mi= ^( 67 r 2 n T ) 2/3 , 

where £ ~ 0.42 is a universal number that was evaluated both experimentally 
EH EH SOI S21 [29J and theoretically [331 13 S3 H2- Integrating Gibbs-Duhem 
identity (see appendix), one then obtains for the imbalanced system 

P) P sr ._L(^) s,2 (f , t + „,)»/*, 

hence P SF /Po = (1 + vf /2 1 '(2£) 3/2 • 

The variation of the pressure versus 77 is displayed in Fig. [TJ We see that for small 
imbalances, i.e. r] smaller than j] c — (2£) 3 / 5 — 1 ~ —0.10, the fully paired superfluid 
is more stable than the fully paired normal phase, confirming the stability of pairing 
against a small mismatch of the Fermi surfaces. The experimental results presented in 
ref. [28] suggest that the two classes of states we have until now restricted ourselves are 
not sufficient to fully capture the physics of imbalanced systems. In particular, a mixed 
normal phase, containing atoms of both species in unequal proportions, must be taken 
into account. A sketch of g(rj) for this intermediate phase is shown in Fig [TJ On this 
more general phase diagram, the parameters r] a and 77/3 are of special importance, since 
they characterize the phase transitions between the three different phases. A glance at 
Fig. [D shows that they must satisfy the inequality 



Va <Vc< Vp, 

and the next section is devoted to an improvement of these bounds. 
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Fig. 1. - Sketch of the grand potential Q as a function of 77 = /^//if. Q is normalized to the 
grand potential of the pure ideal gas of chemical potential /if. Full line: paired superfluid; 
dotted line: fully polarized normal phase; dashed line: intermediate mixed phase. r/ a and 77/3 
designate the critical values for the two transitions between the superfluid/mixed phase and the 
mixed phase/fully polarized Fermi gas. 



3. - The N+l body problem 

Theoretically, the existence of the intermediate phase can be demonstrated by the 
study of the N+l body problem, in other word the study of the ground state of the 
majority Fermi sea in the presence of a single minority atom. This particular system 
corresponds to an intermediate phase with 77 — > r]^ and we will prove that it yields the 
inequality rjp < rj c . 

To address the N+l body problem, we use a variational scheme, that we will compare 
to recent predictions based on Monte-Carlo simulations [46J. Let us consider the following 
trial state 

\i>) = 4> \FS) + J2<t>k, q \k,q), 

k,q 

where |FS) is a spin up Fermi sea plus a spin down impurity with momentum, and 
\k,q) is the perturbed Fermi sea with a spin up atom with momentum q (with q lower 
than hp) excited to momentum k (with k > hp). To satisfy momentum conservation, 
the impurity acquires a momentum q — k. 
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The energy of this state with respect to the non interacting ground state is (H) 
(H ) + (V), with 



(l/>\Ho\l/>) = l^,q| 2 ( 6 fc + e Q~k ~ e q ), 

k,q 

and 



(il>\V\i/>) = y EW 2 + E <t>»«<t>l, q + £ + + ^o«,,) 

\ 9 k,k' ,q k,q,q' q,k 

where e k = h 2 k 2 /2, and the sums on g and k are implicitly limited to q < kp and k > kp- 
As we will check later, <j) k ^ q ~ 1/k 2 for large momenta (see below, eqn. ([TO]) ), in order to 
satisfy the short range behavior 1/r of the pair wave function in real space. This means 
that most of the sums on k diverge for k — > oo. This singular behavior is regularized by 
the renormalization of the coupling constant qb using the Lippman-Schwinger formula. 
It implies that qb vanishes for large cutoff, thus yielding a finite energy. However, it 
must be noted that the third sum in (ip\V\i/j) is convergent and when multiplied by qb 
will give a zero contribution to the final energy and can therefore be omitted in the rest 
of the calculation. 

The minimization of (H) with respect to </>o and (\>k,q is straightforward and yields 
the following set of equations 



(8) yE^ + yE^ = * 

q q,k 

(9) (efc + e q - k - e q )(j) k , q + ^ ^2 4>w,q + ^r</>o = 

k' 

where E is the Lagrange multiplier associated to the normalization of |^), and can also 
be identified with the trial energy. Let us introduce x(q) — + ^Z k 4>k,q- We see from 
eqn. [9] that 



(10) 0k,q 



V e k + e q - k - e q - E 



As expected, we note here the l/e k ~ 1/k 2 dependence for large k. Inserting this 
expression in the definition of x, we obtain 

9b X(Q) 



xM=^o-fE 



V e k + e q - k -e q -E* 
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that is 



/ \ 0o 9b 

x(q) = — 

v 



9B V ^k>k F e fe +e q _ fe -eq 



Finally, eqn. ([H]) can be recast as Efo/gs = ^Z q< k F x(q)/Vi that is, using the explicit 
expression for x(q) : 



E V ^ + 

n<rh^ Or V ' 



q</c F os ^ V Z-^k>k F e h +e q - h -e q -E 

We get rid of the bare coupling constant gs by using the Lippman-Schwinger equation, 
which finally yields the following implicit equation for E 



(11) E=±^ 



V ^ m J_ \^ _J_ 

q<k F 4-Kh 2 c 



V Y^k<k F 2e k + V ^k>k F ( Kek + eq _ k - eq -E 2e k ) 



Before addressing the unitary limit case, let us show that this formula allows us 
to recover the known exact results in the limit of small scattering lengths where the 
denominator is dominated by the 1/a term. The correction to the energy is therefore 



1 47Th 2 a _ 47Th 2 a N 
V ^-f m m V 

q<k F 

where N is the total number of majority atoms. We thus see that the trial state recovers 
the mean-field prediction for low interactions. For a — > + (BEC regime), a little more 
involved calculation allows one to recover the classical molecular binding energy E ~ 
—h 2 /ma 2 . Finally in the case of the unitary regime relevant to experiments, eqn. ([TT]) is 
solved numerically and yields E ~ — 0.3ft 2 /c|i/m, that is rjp < —0.60, a value remarkably 
close to that obtained in Monte-Carlo simulations [46] . 



4. — Trapped system and comparison with MIT experiment 

The model presented in the previous section adresses only the situation of a homo- 
geneous system and to compare with experiments, we need to extend the formalism 
developed in the previous section to the case of trapped systems. To this purpose we 
make use of the Local Density Approximation (LDA), in which we assume that the 
chemical potential of species a varies as 



(12) 



Mr) = M° - V(r), 
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where V is the trapping potential. From this relation, we see that varying r is equivalent 
to varying the chemical potentials of the two species, and in particular their ratio rj(r) 
The two phase transitions described in the previous section will happen for radii r = R a ,/3 
such that iAi(R ay p)/ /j,^(R aj p) = r) a ,(3- Moreover, since the outer rim is constituted by a 
normal ideal gas, the boundary R^ of the majority component is given by the condition 



In an isotropic harmonic trap, we can combine these three relations to eliminate the 
parameters $\ thus obtaining the general relation relating the three radii R a ,f3^- 



where q = (rj a — r/^)/(l — r]p). One striking consequence of this equation is the prediction 
of a threshold at which R a vanishes, corresponding to the disappearance of the fully 
paired superfTuid. This transition happens when the ratio (Rp/R^) 2 reaches the critical 
value q. From the upper and lower bounds obtained for r] a and 77^, we see that q > 0.30. 

This prediction of LDA is remarkably well verified in MIT's experiments [28] for 
which the three phases discussed above were indeed observed, and eq. ([T3]) could be 
tested experimentally (Fig. [2]). On this graph, we see that for large imbalance, the linear 
scaling predicted by eq. ([T3]) is indeed satisfied, with q ~ 0.32, in agreement with the 
lower bound obtained earlier. The deviation from theory observed for (Rp/R^) 2 > 0.5 
is not yet fully understood. However, it must be noted that the discrepancy takes place 
in a regime of low imbalance, where the phase transitions take place in the tail of the 
density distribution. In these regions of low density, we may observe a breakdown of the 
LDA, or of the hydrodynamical expansion that was used to infer the experimental radii. 

The value q ~ 0.32 obtained from the comparison with experimental data can help 
us improve the bounds for rj a ^. Indeed, this relation fixes the relative values of r] a and 
rjp. When combined with the bounds found in the previous section, we obtain indeed 



From the previous analysis, we see that the combination of theoretical arguments and 
analysis of experimental data allows for a precise determination of the thresholds of the 
different phase transitions. Knowing the values of rj a ^ as well as the exact equation of 
state in the fully polarized and fully paired phases, we can even obtain some upper and 
lower bounds for the equation of state of the mixed phase, using the concavity of the 
grand potential [36] . 



/i T (# T ) = 0. 



(13) 




(14) 
(15) 



0.62 < rip < -0.60 
0.10 <f] a < -0.088 



Unitary polarized Fermi gases 




(Rg/R t )' 



Fig. 2. - Comparison of equation (|13[) with experimental data from MIT. A fit to the data yields 
q ~ 0.32. Inset: sketch of the density profile. The full (resp. dashed) line corresponds to the 
density of the majority (resp. minority) component. R a marks the end of the superfluid region, 
Rf3 that of the mixture and is the frontier of the majority cloud. 



5. — Elongated systems and Rice's experiment 

Surprisingly, similar experiments performed at Rice University showed no evidence of 
an intermediate phase, but rather the coexistence of the fully paired and fully polarized 
phases only. Measurements of the axial radii of the two phases from ref. [29] are presented 
in Fig. [3] and can be compared with the model presented above when omitting the 
intermediate mixed phase [37]. In these conditions, the inner superfluid region is now 
defined by the condition /i|(r)//i|(r) < r] c and is bounded by the radius R± denned by 



,16) fl;= 4,(4^Y 

Atoms of the minority species are located in the paired superfluid phase only. We 
thus have 



(17) 



J 

J r<R\ 
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Fig. 3. - Rice's radius measurement and comparison with a two phase model. The radius Ri is 
scaled in units of the Thomas Fermi radius of an ideal gas with a the same atom number Ni. 



where R 2 = (fi® + [i®)/muj 2 and 
(18) 



x \J\-x 2 (-3 + 14 x 2 - 8 x A ) + 3 arcsinO) 
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Excess atoms of the majority species are located between r = R± and r = R^ such 
that mO B%/2 = The number of excess atoms is therefore iVj — = J R T n\{r) d 3 r, 
hence 



(19) 



^T-^l = |^fe) (ff(l)- 5 (i?i/ii T )). 



Dividing by ([T9|) by (fTTj) yields the implicit equation for rjo = /jP^/ M| as a function of 



(20) 



^T_ 1 ,,3/2 8 g(l)- g (^/iZ T ) 
Ni ? (l + r?o) 3 s(i?i/i?) ' 
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Equation (|20|) is solved numerically and the value obtained for 770 is then used to 
calculate the radii and R± . The predicted evolution of the Ri versus the population 
imbalance P = (7V T - iV|)/(iV T + N±) is shown in Fig. H To follow Ref. [29], we 
have normalized each Ri to the Thomas-Fermi radius Rtf associated to an ideal gas 
containing Ni atoms. The agreement with the experimental data is quite good as soon 
as P > 0.1, a remarkable result, since the model presented here contains no adjustable 
parameter, as soon as the value of £ is known. 

Despite this remarkable agreement, this simple two phase+local density approxima- 
tion model fails to captures all experimental features. In particular, a qualitative discrep- 
ancy occurs in the comparison between the theoretical and integrated density profiles. 
Indeed, as shown in [37]? LDA at unitarity implies a constant density difference in the 
paired superfluid region, in contradiction with experimental data. One solution to this 
problem was presented in [32[ [31]. In these papers, it is noted that in the presence of 
phase transitions, the description of the sharp frontier separating to adjacent phases in- 
volves the introduction of density gradient terms in the energy. When the interface in 
thin enough, they can be encapsulated in a new surface tension energy term reading [3T] 



Js 

where S is the interface between the two phases, and 7 is the surface tension constant, 
which should dimensionally vary as 



Here, A is a numerical factor that will be determined by comparison with experi- 
ments and we have used the fact that at the coincidence between the phases, the ra- 
tio is fixed and equal to rj c , meaning that the two chemical potentials are no 
longer independent. We can minimize the total grand potential ft = Qbuik + ^st, where 
^buik = — J (Pn + Psf) d 3 r is the bulk contribution to the energy. Following [31] we 
simplify the analysis by assuming that the interfaces are ellipsoidal, and for A ~ 10 -4 
one obtains the results presented in Fig. 0J which coincides with experimental data. The 
absence of capillary effects at MIT can be explained by a smaller trap aspect ratio and 
a larger atom number of atoms compared with Rice's experimental situation, as shown 
by a simple scaling argument [31]. 

6. — Conclusion 

The formalism presented here allows for a simple and quantitative description of 
macroscopic properties of polarized Fermi gases in the regime of strong interaction. This 
analysis is nevertheless far from being complete, since it does not give any information on 
the superfluid nature of the various phases. For instance, the mixed region of the phase 




7 = A 





Fig. 4. - Integrated density difference in Rice experiments, and comparison with the surface 
tension model (data from [31]). Dashed line, LDA prediction: the density difference is fiat in 
the superfluid region, in contradiction with experimental date. Full line: Two phase model 
incorporating surface tension effect. The same surface tension parameter A =~ 10 -4 is used for 
all three graphs. 



diagram may contain superfluid and normal subdomains, the transition between this 
two regimes being characterized by a universal number r? 7 e [r]p,r) a ]. The quantitative 
understanding of these superfluid properties will require beyond mean-field theories, such 
as the Monte-Carlo calculations of [46] . 
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8. — Appendix: t her mo dynamical relations for the grand potential 

Let us consider a homogeneous many-body system characterized by a hamiltonian Ho 
and containing particles of p different species labelled by i = l..p. In the grand canonical 
ensemble, one looks for the ground state of this system by letting the atom numbers 
fluctuate, but keeping the expectation values (Ni=i, mP ) fixed. This therefore requires to 
find the ground state of the grand potential = H — J2i=i MiM, where the Hi are 

Lagrange multiplier that we interpret as chemical potentials. 

Let \ip(fjLi,V)) be the ground state of the grand potential, we set Q (/i^ , V ) = (ip (fii , V ) | O \ip (/i^ , V ) ) . 
Using Hellman-Feynman relation, we can write that 



(21) |^ = MlH,V)Ail,OH,V)) = -(Ni. 

from which we deduce that 
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(22) dQ = -Nidfn + —dV. 

By definition, and by analogie with classical thermodynamics, we identify dyQ with 
— P, the pressure in the system. 

Let us now us the extensivity of the potential: when the volume is multiplied by 
some scaling factor A, Vt is multiplied by the same factor. In other words, we have 
ft(AV,/ii) = Aft(V,/ii). Taking A = 1/V, we get ft(V,/ii) = Differentiating 
this with respect to V, we note that f2(l,/Zi) = — P, hence 



(23) n = -py 

From this equation, we see that the minimum grand potential is state has also the 
highest pressure. P can moreover be calculated by differentiating Vt and using equations 
([23]) and ([22]) . We then obtain the Gibbs-Duhem relation 



(24) dP = J2 n i d ^ 

i 

where rii = Ni/V is the density of species i. From equation ([24]). we see that the pressure 
(hence the grand potential) can be obtained simply from the knowledge of the equation 
of state rii(fij). 

8*1. Concavity. - Since, by definition, IVKaO) * s the ground state of tlfai), we have 
for any 5/iii 

(25) (V(/i; + ^)l^(^)l^(^ + fyi)} > 

Moreover, if one notes that fl(fii) = + + 5^ ■ SfijNj, we see that for any 

(26) + + y ^5fijN j (iii + J/ii) > 

Finally, recalling that = <9 M . ^ and after expansion of equation ([26]) to second order 
in SfjLi, we obtain 



hence proving the concavity of the grand-potential (or conversely the convexity of the 
pressure). 
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